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Abstract 

The rupture of rubber differs from conventional fracture. It is supersonic, and the speed is 
determined by strain levels ahead of the tip rather than total strain energy as for ordinary 
cracks. Dissipation plays a very important role in allowing the propagation of ruptures, and 
the back edges of ruptures must toughen as they contract, or the rupture is unstable. This ar- 
ticle presents several levels of theoretical description of this phenomenon: first, a numerical 
procedure capable of incorporating large extensions, dynamics, and bond rupture; second, 
a simple continuum model that can be solved analytically, and which reproduces several 
features of elementary shock physics; and third, an analytically solvable discrete model 
that accurately reproduces numerical and experimental results, and explains the scaling 
laws that underly this new failure mode. Predictions for rupture speed compare well with 
experiment. 
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1 Introduction 



The theory of fracture was originally developed to explai n the failure of brittle ma- 



terial s , where cra cks have a number of common features (Irwin, 1957; Kanninen and Popelar 



1985 : Thomson , 1986h . A stress singularity builds up in the vicinity of the tip. 



Stresses diverge as 1/y/r, where r is the distance to the tip. The displacement of 
material near the tip behaves as y 7 ?, which means that the tip viewed closely has the 
shape of a sideways parabola. Energy flows in towards the crack from far away and 
concentrates itself at the tip in just the amount needed to snap bonds and feed other 
dissipative processes. Partly for this reason, cracks in tension cannot travel faster 
than the Rayleigh wa ve speed, which is t he speed at which sound travels across a 
free surface(Broberg, 19991: Freund . 1990h . 
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Experimental evidence has accumulated showing that the rapid rupture of rubber 
sheets, such as when one pops a balloon, is different. If one cuts a horizontal slit 
in a rubber sheet and stretches it up, the opening profile is a sideways parabola, as 
expected for static cracks. Once the rupture begins to run, however, the character- 
istic \fr opening dis placement disappears, and is replaced by a wedge-like shape 



( Deegan et all 120021). Furthermore, the ru pture speed exceeds the shear wave speed 



upt 
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in advance of the tip (Petersa n et alll2004l) . However the extensions in rubber when 



it ruptures are very large. The ordinary theory of fracture begins with the assump- 
tion that strains are very small, typically on the order of a fraction of a percent 
far ahead of the crack. In rubber, rupture initiates when strains are on the order of 
several hundred percent. Therefore it is not clear how much of fracture mechanics 
ought to apply to rubber, and whether the violations of rules about rupture speed 
are simply the natural result for a material that breaks for very large extensions, or 
whether the mode of failure is something new. 



In this article I will explain the case first outlined in Marder (2005) that the rupture 
of rubber is different from conventional fracture. It is a tensile failure. The ruptures 
always travel faster than the shear wave speed. The opening is a wedge that obeys a 
simple relation that applies to Mach cones. Stress is singular near the tip, but strain 
is not, and material in front of the tip must be brought rather near the point of failure 
for ruptures to propagate. 

Buehler, Gao, and Abraham (2003) have proposed that hyperelasticity plays a crit- 
ical role in dynamic fracture, and that in particular an increase of sound speed near 
a crack tip can allow cracks to travel faster than the distant shear wave speed. In 
the analysis of this paper, there is an increase of elastic modulus near the tip of the 
rupture, but it is not in the form that Buehler, Gao, and Abraham proposed. The 
theory here relies upon an increase in the modulus of rubber with frequency, rather 
than upon the increase in modulus of rubber with extension. The low-frequency 
modul us of real rubb er certainly increases as rubber stretches towards the breaking 



point (Treloar, 1975, p. 2). However, I found in numerical investigations that the 
behavior of ruptures does not change appreciably whether such stiffening prior to 
breakage is included or not. In the computations of this paper, Lagrangean sound 
speeds are either independent of extension (Sections 5 and 6), or else increase as 
the rubber contracts (Section 4.3). In view of the fairly detailed comparison of the- 
ory, numerical work, and experiment obtained in this work, the conclusion is that 
static hyperelasticity is not relevant to the supersonic rupture of rubber. 



The theory comes in several forms. I begin with a numerical model of rubber that 
is fairly realistic and includes most of the physical features of rubber indicated by 
experiment. There are only three parameters in this model not determined directly 
from experiment, and of those there is only one to which the fracture dynamics are 
particularly sensitive. Next, I note that after stripping some of the realistic com- 
plexity out of the numerical model, the dynamics of ruptures scarcely change. The 
resulting theory is so simple that it can be solved analytically. The analytical solu- 
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tion takes two forms. The first is a continuum model with a simple failure criterion 
that leads to compact closed-form expressions for rupture speed and opening angle. 
The expressions for rupture velocity agree with laboratory data within experimen- 
tal error, although they disagree with rupture speeds from numerical modeling by 
around 10%. Finally, the discrete rupture theory has a complete analytical solution. 
This solution enables one to see the relation between the conventional theory of 
fracture and supersonic ruptures, and is particular how these two types of solution 
scale in the macroscopic limit. 

The structure of the paper is as follows: Section 2 presents some elementary physi- 
cal ideas to describe rubber rupture. Section 3 lays out the continuum energy func- 
tional for rubber on which much of the subsequent discussion will be based. Section 
4 develops a computational method that makes it possible to obtain numerical solu- 
tions that mimic the main features of experimental ruptures. Section 5 extracts from 
the numerical model a simple continuum theory and presents its solution. Section 6 
proceeds further to obtain an exact analytical solution of the full numerical model, 
after some simplifications. Section 7 compares theoretical predictions with the ex- 
perimental results. There are also five appendices, which discuss (A) how to obtain 
an effective two-dimensional theory from the original three-dimensional theory, 
(B) the computation of sound speeds from the two-dimensional continuum the- 
ory, (C) the computation of forces for the numerical model, (D) a lattice instability 
found in some of the numerical models, and (E) the solution of the discrete model 
by Wiener- Hopf techniques. 



2 Elementary considerations 

Much of this paper will be concerned with the speed of ruptures in rubber. There- 
fore it is necessary to begin by describing precisely how speed will be defined. 
Suppose that a sound wave or a rupture travels through a highly stretched rubber 
sheet. One can choose either to describe its speed in the laboratory (Eulerian sound 
speed) or back in a coordinate system tied to the original location of mass points 
(Lagrangean sound speed). To be more explicit, consider a sheet of rubber lying 
relaxed in some initial material reference configuration. Describe the sheet in this 
configuration with the variable f, which serves as a label for mass points. Once the 
rubber is stretched and begins to move, the new locations of mass points will be 
given by the variable u. Now consider some deformation or bump moving through 
the rubber at constant speed. The maximum amplitude of the bump is located at 
some point u(t) which changes in time, and the speed at which it moves is the lab- 
oratory velocity. However, there is another velocity, which from a theoretical point 
of view is much more natural to employ, and which helps assemble the experimen- 
tal results into a compact scaling form. Whenever the bump is located at u, one can 
identify the original location f of the mass point now at the center of the bump, and 

— * — * 

r(t) also evolves in time. The speed at which r(t) travels is the Lagrangean speed, 
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Figure 1 . Illustration of material and laboratory frames used to describe rupture of rubber. 
The top left panel shows a rubber sheet before it has been stretched, and the top right shows 
it stretched by amounts X x and X y in the x and y directions. The lower right panel shows a 
traveling rupture as seen in the laboratory, while the lower left shows the same rupture back 
in the material frame. If the normal velocity of the rupture is c, then the forward velocity 
of the rupture must be v sin 9. The shock lines shown in the lower left panel correspond to 
the leading edge of compressed material in the lower right panel. 

and unless otherwise specified, sound and rupture speeds will always be measured 
in this Lagrangean reference system. For example, if one considers a bump of small 
amplitude moving along x in rubber that has been stretched by a factor of \ x above 
its original length, then the speed of the bump in the laboratory frame exceeds that 
in the reference frame by a factor of X x . 

Now consider a thin sheet of rubber that is stretched by factors of \ x and \ y in 
the x and y directions respectively, as shown in Figure 1 . Stick a pin into the sheet 
on the left hand si de so that a rupture runs to the right along the x direction. One 
sees in experiment ( Petersan et al. . 2004 : Deegan et al.l E002) that the rupture con- 
sists in two straight fronts that meet at a point, forming a wedge. Suppose that the 
two straight fronts are shock fronts, traveling at th e Lagrangean wave speed c. As 
is customary for the elementary theory of shocks ( Serway an d Beichnet . 200d p. 
534), the speed v of the tip of the rupture must obey 



c 
v 



sin.6 1 , 



(1) 



where 9 is the opening angle of the rupture, as shown in Figure 1 . The slope of the 
upper face of the shock line is — tan 9, which is 



material frame slope 



^Jv 2 /c 2 - 1 ' 



(2) 



In the laboratory, where distances are stretched by factors of and X y , the slope 
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will intead be 

laboratory slope = . = =, (3) 

A.yV/c 2 - 1 

since if one draws a line of slope a on a sheet of rubber and stretches it along x by 
a factor of X x , the slope decreases by a factor of A^. 

Eq. (3) will re-emerge from detailed calculations as Eq. (69). The main physical 
quantity left undetermined is the rupture velocity v. A decent approximate relation, 
obtained as Eq. (79) that closes the theory is 

\} = \\l + \\l/{l-c*/v% 

where A/ is an extension at which polymers in rubber snap. Far ahead of the rupture, 
the bonds that will eventually be brought to the snapping point are already stretched 
an amount y^A 2 + |A 2 over their original length. In simple physical terms, then, 
the assertion is that just in front of the tip of the rupture, material stretches by an 
additional factor of 1/(1 — c 2 /v 2 ) in the vertical direction. However, I have not 
found an elementary argument to produce this relation. 



3 Continuum Energy Functional 



3. 1 Coordinate system and definition of energy functional 



Strains in rubber are several hundred per cent at rupture and one must use n onlinear 
elastic theory to describe the situation (Atkin and Fox . 1980: lQgdenL [l9R4l).I state 
just enough of the theory to establish notation. Adopt a description of a highly 
deformed rubber sheet with 



\ r xi r y) 



u 



(U , u 



(4) 



The original location of all mass points is given by r and the location of points after 
the rubber is moved and stretched is given by u. Note that u is measured from the 
origin, not from the original location of the mass point r. Define the Lagrangean 
strain tensor 



E a /3 = \ 



^ du 1 du 1 



(5) 



., dr a dr p 

From this strain tensor one can define three rotationally invariant quantities. These 
are 
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lf D = Tr E 



r 3D 



a</3 



\E aa Ej3i3 — E\ 



/f D = det£, 



(6a) 
(6b) 

(6c) 



Rubber is highly incompressible JTre1oaiirr975[ p. 61). Accordingly, for a thin sheet 
of rubber, one can express the thickness at every point in terms of the strains in the 
x — y plane, and project the theory into two dimensions, as discussed in Appendix 
A. In two dimensions one has only the two invariants, 



Ii = Tr E; I 2 = E xx E yy — E 2 xy , 

and using the incompressibility of rubber to solve for E zz one finds 

1 



E z 



4/ 2 + 2/i + 1 



1 



(7) 



(8) 



in terms of which the first two of the three-dimensional strain invariants take the 
form 



I\ D = h + E zz 
1 1° — h + E zz Ii 



(9) 



The energy density of a thin sheet of rubber is then taken to be given by some 
function e of Ii and 7 2 , and the total energy U is 



U = p d? e(h{r'),I 2 {r')) . 



(10) 



The integral is performed in the material frame, and the mass density p is also 
measured in the material frame. 

An energy-conserving equation of motion for this theory is 



5 



5ttT(r) 



dr 1 \p 



u 



pe 



(11) 



where p is the mass per area, again measured in the material frame. Performing the 
functional derivatives, one has 



pu 1 



_5U_ 

8ift{r) 



(12) 



Appendix A demonstrates that this equation of motion for a two-dimensional sheet 
obtained by calculating E zz through Eq. (8) is the same one obtains from the 
Piola-Kirkhoff stress tensor after imposing incompressibility and requiring that the 
Cauchy stress T zz vanish. 
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3.2 Sound Speeds 



The experiments by Pet ersan et al.l (120041) that stimulated this study obtained de 



tailed information about the speed of sound in rubber under a range of loading 
conditions. For a while, we found the results puzzling, but eventually realized that 
they could all easily be explained by the Mooney-Rivlin theory. Appendix B con- 
tains a sketch of how to obtain sound speeds from the equ ation of motion (12). 
Here I record only the final results, all of which are standard(Eringen an d Suhubl 



1974, v. 1, pp. 120, 263). Suppose that a rubber sheet is strained uniformly with 



displacement field 

u = (X x x + s xy y, X y y + s yx x). (13) 

Look for the speed of sound along the x and y axes of a sample that is extended by 
the two factors X x and X y ; s xy and s yx are included in Eq. (13) only because one 
must be able to perform calculations involving small virtual shears around this base 
state. Then there is a longitudinal sound wave along the x axis whose speed with 

$xy Syx is 



2 _ 

° xl ~ d\l' 

Similarly, the speed of longitudinal waves in the y direction is 



(14a) 



d 2 e 

c « - M (14b) 

There is also a shear wave that travels along x and is polarized along y with speed 



4 = |f (wo 

us yx 

Similarly a wave traveling along y and polarized along x has speed 



4s = q^-- ( 14d ) 



Alternatively, one can express sounds speeds in terms of derivatives with respect to 
strain tensor components. One has for the longitudinal wave speed along x, 



2 de 2 d 2 e 

° lx ~ 8E~ + x 'dE r ( } 



1 



while for the shear wave speed (setting E yx = E xy before taking the derivatives) 



yx — ^xy 

c 2 = ^3 — (15b) 
dE xx + AdEl y U3D; 



All of these speeds are Lagrangean speeds, as described in Section 3.1. 

Sound speeds provide a convenient way to assemble experimental data about the 
constitutive behavior of rubber. In some cases, sound speeds are measured directly 
through time of flight, while in other cases they are measured through small ex- 
te nsions of the sample a round a base state as suggested by Eqs. (14). The results 
of IPetersan et alJ (120041) are quite simple. Over a range of biaxial states where 



\ x £ [2, 3.5] and \ y £ [2, 3.5] the Lagrangean wave speeds appear to be constant, 
with the longitudinal wave speed around 20% greater than the shear wave speed. 
From Eqs. (15) one finds that this is not possible. The only way for the longitudinal 
and shear wave speeds both to be constant is to adopt e(Ji, I2) oc h, and in this 
case the longitudinal and shear wave speeds must be equal. 



This apparent difficulty is resolved by exam ining a bit more carefully the free en- 
ergy fu nctional due to Mooney and Rivlin dTreloai . 1975 : Moonev . 1940l: RivlinL 
194833)- 

The Mooney-Rivlin theory says that the free energy density of rubber is 



U/p 



3D 



Bl 



3D\ 



(16) 



where U has units of energy per volume, p is mass density, A is a constant with 
units of velocity squared, and B is dimensionless. Using Eq. (9) one obtains an 
effective two-dimensional Mooney-Rivlin theory 



e{h, h) = A(h + BI 2 + E zz (l + BI 2 )) 



(17) 



For extensions \ x and X y on the order of 2 or greater, E zz + 1/2 is of order 1 / (\ x \ 2 ^ 
and is at least 64 times smaller than E xx or E yy . Therefore, for the purpose of 
examining the experiments, it is sufficient to use 

e(Ji, h) = A(h + BI 2 ) = A [(E xx + E yy ) + B [E xx E yy - E 2 xy )] (18) 

Employing Eqs. (15) and (18)one finds for longitudinal and shear wave speeds 
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(19a) 
(19b) 
(19c) 
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Figure 2. (A) Sound speed data for various values of A^, and X y = 3.2. Experimental 
longitudinal (c x i = □) and shear (c xs = ■) speeds are shown in the material frame. They 

are roughly constant, and can be fit by (19a) ( — ) and (19c) ( ), using A = 501 (m/s) 2 , 

B = 0.106. The shear wave speed comes out to 21.8 m/s and the longitudinal wave speed 
for X y = 3.2 is 27.3. (B) Using the constants A and B obtained from the data in (A), 
calculate longitudinal wave speeds in the y direction c y i. According to Eq. 19, this is the 
only speed that varies as a function of \ x . The agreement with the data is satisfactory. 

Thus for a Mooney-Rivlin material the shear wave speed is constant, and the longi- 
tudinal speed along x is independent of X x but depends quadratically on X y . Turning 
to the experimental data, one finds that they are consistent with these observations 
as shown in Fig. 3.2, and that one can fix the constants A and B. 



For studying the rupture of rubber, the energy density in Eq. (17) is both too simple 
and too complicated. It is too simple because it does not account for the fact that 
when rubber is stretched enough, the polymers pull apart and the force between 
adjacent regions drops irreversibly to zero. It is too complicated because the terms 
involving I 2 and E zz produce nonlinear equations of motion that are impossible 
to solve analytically. Therefore, to analyze the problem, I will pursue two differ- 
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ent routes. First, I will discuss numerical routines that supplement Eq. (17) with 
information about rupture, toughening, and dissipation, and produce supersonic 
solutions. Second, I will isolate from Eq. (17) terms that are sufficient to produce 
good agreement with numerics and experiment, while simplifying matters enough 
to permit analytical solution. 



4 Numerical methodology 



For the problem of fracture, there are great advantages to thinking in terms of 
molecular dynamics. From a continuum viewpoint, it is difficult to understand how 
to construct a physically sensible theory where material gives way. From an atomic 
viewpoint it is easy; when two atoms are separated by more than a certain distance, 
they stop applying force to one another. Therefore, I have found a simple set of 
microscopic interactions that produces the Mooney-Rivlin theory of Eq. (17) in 
the continuum limit. The interacting mass-points that appear in the theory should 
not be thought of as atoms. To describe rubber, they should be thought of as nodes 
in a cross-linked polymer network, with a characteristic spacing of around a mi- 
cron. There are some possible objections to this approach. Rubber is much more 
complex than a triangular lattice, mass is distributed rather than being concentrated 
at nodes, and one might worry about the fact that the two-dimensional array of 
nodes has been engineered to reproduce dynamics that derive from projections into 
two dimensions of three-dimensional equations of motion. There is no complete 
answer to these objections; the best response is to show that the resulting theory 
provides detailed correspondence with experiment, and allows much analytical and 
numerical progress. 

Simi lar numerical techniques have been used before; for example by Seung an d Nelson! 
(1988). The techniques of that paper must be extended to include bond snapping 
and dissipation, which I carr y out here. The philo s ophy is also similar to the Vir- 
tual Internal Bond method dGao and Kleinl 1 19981: iKlein and Gaol [l998) and the 



Peridynamic Model dSilling and BobaruL 120051: ISillingL l20Q0h . which also focus on 
a collection of discrete interacting mass points considerably larger than atoms in 
order to obtain rules for fracture. However, all details of the implementation of 
this idea are different; the formulation presented here has the advantage of leading 
in one case to a discrete model of nonlinear materials with a complete analytical 
solution. 



4. 1 Low-order polynomial terms for microscopic theory 



Consider the triangular lattice depicted in Fig. 3. The figure shows the original 
locations of all particles prior to any distortion, denoted by a i} and the lattice spac- 
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Figure 3. Diagram showing triangular lattice of lattice spacing a and nearest-neighbor 
vectors dlj used in this section. The original locations of particles at rest are given by a, 
while their location in the laboratory after deformation is given by u. 

ing is a. After distortion, the position of particles in the laboratory is given by 
UjThe goal is to construct a theory for the energy required to displace particles on 
this lattice that involves I x and J 2 , and employs quadratic and quartic functions of 
displacements u. It is possible to construct such a theory by considering simple 
combinations of rotationally invariant operations on nearest-neighbor vectors. Let 
let nil) refer to the nearest neighbors of i, and define 



/ - >' {•;::' : " g2) ifu ^ <x f ( 20a) 

else 



1 x . Uuij ■ Uij - a 2 ) 2 ifuij<X f 




jen(.-)l(A/-a 2 ) else 



Hi = E H u ij) h ( u ik) (uij ■ u ik + 2a 2 ) 2 , (20c) 
and %) = 1/(1 + e^ A ^ Us )- (20d) 



The sums are carried out over the 6 nearest neighbors of point i shown in Fig. 3. 
The terms are constructed so as to become constant and therefore describe breaking 
bonds when increases to more than A/. The final term requires a cutoff function 
h since this is the only way to ensure both that Hi be continuous, and that it settle 
down to a constant value when or are large; the constant u s describes the 
scale ov er which h vanishes. Terms of this form are standard in molecular dynam- 



ics (e.g. Stillinger and Weber (1985)). To form a correspondence with continuum 



theory, suppose that no bonds are stretched past the breaking point (Uij < A/), and 
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approximate the position of neighbors of point % by 



Insert Eq. (21) into Eqs. (20a) to obtain (using Eq. (5)) 



1 / q du a du a „ 

- > > aj,-— — — — a , — a 

jen(t) \op7 p i 



\ E E<[2% + W«S-« 2 

je(i) \ /37 
(■E'xx + E yy )a 2 = ha 2 , 



- 9 E E« 
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= g E E 4 [2^7 + <W 

jen(i) \a/37 
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du a du a 
' dxp dx 7 



a] k + 2a 



20 

+ -Eyj/) 2 + — - E XX Eyy^j + 4 

J? - + 4l a 4 . 



Comparing Eqs. (24), (22), and (25) with Eq. (7) gives 



F 

Tl — 1 
11 ~ n2 



12 



(28) 



or alternatively, 



9 1 



8a 



4 (G t -H t + A) 



(29) 



Numerically, Eq. (28) is much less costly to compute than Eq. (29). However, Eq. 
(28) has the unfortunate property that under biaxial strain, spatially uniform states 
are unstable when this representation of I2 is employed. Particles bunch up in a 
non-uniform way within each unit cell, forming stripes on a microscopic scale, as 
shown in Appendix C. It could be that this behav ior is related to the physical phe- 



nomenon of strain crystallization (Treloait ll975L p. 20). However, as strain crys 



tallization does not occur experimentally in the range of extensions where ruptures 
are observed, I have largely employed Eq. (29) in preference to Eq. (28). 



4.2 Specification of numerical energy functional 



To form a numerical representation of Eq. (10), take 



U = mJ2e(ll,P 2 ), (30) 

i 

where m is the mass in a unit cell. Then if f2 is the volume of a unit cell, 

u ~!il dre {W),W)), (3i) 

so since m/VL = p, e in Eq. (30)corresponds to e in the continuum theory, and has 
units of velocity squared. In particular, for the Mooney-Rivlin theory, one has 



e(ll P 2 ) = A{I[ + BP 2 + El z {\ + BP 2 ), (32) 

where I{ is given by Eq. (27), I 2 is given either by Eqs. (28) or (29) , and E\ z is 
given by(8), with I\ and I 2 substituted for Ii and I 2 . 



4.3 Equation of Motion 



Given the energy functional (30) one can obtain the force on every particle and 
therefore an equation of motion. In addition to the conservative force resulting from 
derivatives of the energy, add Kelvin dissipation, so that the complete equation of 
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Figure 4. Four panels showing initiation of rupture in numerical sheet of rubber. The nu- 
merical system contains around 70,000 particles, and solves Eq. (33), using Eqs. (30) and 
(32) with A = 501m 2 /s 2 , B = .106, A/ = 5.5, and (3 = 3. The first panel shows the initial 
pop, the second shows the system 12.5 time units later, the third after 25 time units, and the 
final panel after 250 time units, where the time unit is a/\fA . 

motion reads 

muf = -dU/du? + £ *?£lv°e(\ f - Uij ). (33) 

jen(i) 

The final term in Eq. (33) represents the Kelvin dissipation, and is the simplest 
dissipative term permitted by symmetry. The motion of each mass point dissipates 
some energy in proportion to its velocity relative to each neighbor. This dissipation 
vanishes when the bond between two neighbors breaks. 

The computation of dU/duf is not particularly difficult as force computations go. 
A formalism that makes it easy to exploit symmetry to reduce the amount of com- 
putation is briefly described in Appendix C. 

One final rule is employed in the numerical runs, although it is not indicated ex- 
plicitly in Eq. (33). Whenever some bond drops to a length less than 1.5a, the 
failure extension A/ for the remaining bonds attached to nodes % and j increases. 
Without some rule of this type, the back faces of the crack disintegrate. The reason 
for this rule will be explained in Section 5.3. 

Figure 4 shows characteristic panels from a numerical solution of Eq. (33). First, 
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Figure 5. Comparison of experimental crack speeds and speeds in simulation. Simula- 
tions and experiments are conducted with \ x = 2.2, and a range of vertical extensions 
X y .The parameters describing the properties of the continuum nonlinear elastic theory are 
A = 501 (m/s) 2 , B = 0.106, as in Figure 3.2. Bonds in the simulation fail when extended 
an amount A/ = 5.5 above their original length. The magnitude of Kelvin dissipation is 
(3 = 3. In addition, the figure displays results from a substantially simplified numerical 
model where B = 0, and where E zz is set to zero as well. Stripping most of the complexity 
from the numerics has little effect on the results. 

one prepares a uniformly strained sheet, in this case with extensions A^ = 2.2, 
\ y = 3.2. Next, two rows of particles are selected near the left hand side of the 
sample: the rows are five particles wide, and they sit right on top of one another. 
The upper particles are given a large upward velocity, and the lower particles are 
given a large downward velocity. This initial condition has the effect of popping a 
hole in the strip. As shown in the subsequent panels of Figure 4, the hole initially 
develops in a circular fashion, but as it senses the upper and lower boundaries, it 
begins to run sideways, and eventually turns into a shock-like rupture front that 
travels in steady state indefinitely to the right. 

Figure 5 shows a comparison of experimental rupture velocities once steady state 
has been reached with results from numerical simulations. The agreement is satis- 
factory. However, the figure also contains the results of a different set of simula- 
tions. In these, Eq. (33) is solved not with the Mooney-Rivlin energy function in 
Eq. (32), but with a much simplified Neo-Hookean energy functional 

e NH (Iir 2 ) = c 2 I\. (34) 

The velocity of ruptures described by this very simple theory is indistinguishable 
from the velocity of ruptures described by the more elaborate Mooney-Rivlin the- 
ory. This observation opens the way to accurate analytical descriptions of the rup- 
ture of rubber, both at the continuum and discrete levels. 
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5 Continuum Neo-Hookean theory 



5. 1 Cracks in the Neo-Hookean theory 



The continuum Neo-Hookean energy is 

eNH = Alt 



A 
~2 



du x \ 2 ( du x ^ 



dx 



dy 



+ 



dx 



+ 



'du* ' 



(35) 



This energy functional was first employed for the stu dy of rubber by Mooney 
(1940), and was employed for the study of fractures bv lKlingbeil and Shield! d 19661) . 
It has a number of interesting properties. Because it is quadratic in displacements, 
it leads to a linear equation of motion for u, which is easy to approach analytically. 
Nevertheless, it describes very large displacements of rubber, and in that sense is 
still a nonlinear theory. According to Eq. (19), there is only one sound speed for 
this theory, 

c 2 = A. (36) 
The equation of motion of the Neo-Hookean theory follows from Eqs (22), (30), 
and (33), and is 



mil? 



E 



2mc 2 



(37) 



In the continuum limit, overlooking bond rupture, one has 



mm 



E K - <) + W - <) 



2mc 
3a 2 

mc 2 d 2 (u a + (3u c 

a«7 



dr a dr~ 



a 7 
~ a ij a ij 



(38) 



mc 2 V 2 {u a + (3ii a ) 



For the study of fracture, there are some additional simplifications that arise in 
simple geometries. Consider a semi-infinite crack moving along the center line of 
an infinitely long strip as illustrated in Figure 1 . The bottom of the strip is held at 
height — b\ y , while the top of the strip is raised rigidly to height \ y b, where 2b is 
the original height of the strip. Suppose that the rubber is initially stretched by a 
factor \ x everywhere in the horizontal direction. Dropping for the moment Kelvin 
dissipation, the equation of motion is 



u 



c z V z u 



2„.£ 



(39a) 
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yV = C 2 VV 



(39b) 



and the boundary conditions for a crack with tip at vt are 



u y (x, -b) = -X y b; u y (x, b) = X y b; (40a) 
du y /dy\ y=0 = for x < vt (40b) 
u y {x,0)=0 forx>vt (40c) 

The main point to make here is that the boundary conditions nowhere involve 
u x , nor do the equations of motion couple u x and u y . Therefore, one can take 
u x (x, y) = \ x x for all time, and the problem reduces to one involving only u y . 
This problem is mathematically identical to the problem o f crack motion in anti- 
plane shear, which is discussed in textbooks. For example, (Brober 3- ll999L p. 127) 
provides the solution of a stationary crack in this geometry, and the solution for 
a crack moving at steady velocity v can be obtained from the static solution by a 

simple change of variables x — x/ ^1 — v 2 /c 2 . The solutions become increasingly 
blunt as they approach the sound speed c. 

I have not included Kelvin dissipation in Eq. (39). The reason is that this term 
would destroy the conventional y/r singularity expected for cracks. If one supposes 
there exists a solution where u y (r, ) ~ y/r, then Kelvin diss ipation produc es an 
infinite amount of energy dissipation in the vicinity of the tip(Marder, 2004). We 
will see, however, that for supersonic solutions where the mathematical structure 
near the tip is different, Kelvin dissipation is not only permitted but required. 



5.2 Shocks in Neo-Hookean material 



I now proceed to study supersonic solutions of the Neo-Hookean theory in the 
presence of dissipation. Adopting once again the geometry of Figure 1, one can 
conclude that u x (x, y) = X x x at all times, and focus only upon u y . Since this is the 
only variable to consider in the following discussion, u will refer to u y . The vertical 
displacement u = u y obeys the equation of motion 



u = c 2 X7 2 u + c 2 /3V V (41) 

with boundary conditions at the top and bottom of the strip still given by Eq. (40a), 
but now at y = 

du d 2 u 

— = -(3—— for x < 0; u = for x > 0. (42) 

ay otoy 
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The boundary condition is obtained heuristically by discretizing the derivatives in 
the y direction, and eliminating the near- neighbor interactions for y < 0. That is, 
write 

j u(x + a, y + u(x — a, y) 
V 2 w^ +M (a; ?2/ + a ) +u ( x ^ y - a y (43) 

-Au(x,y)]/a 2 

On the boundary there are no particles located at y — a, then one has there instead 



[ u(x + a, y + u(x - a, y) - 2u(x, y) 
V u ks . (44) 

+ {u(x, y + a) — u(x, y)}]/a 2 

The term in curly brackets must vanish, or it produces contributions of order 1 /a. 
Analyzing also the last term of Eq. (41) in this way produces Eq. (42). 

In steady state, the equation of motion and boundary condition become 



^ 2 £ = c 2 V 2 M -^c 2 V 2 ^, (45) 
with boundary condition at y = 

du d 2 u 

— = v(3 - - for x < 0; u — for x > 0. (46) 
ay oxoy 

This system can be solved by the Wiener-Hopf technique. Consider Eq. (45) for 
y > 0. Subtract out the asymptotic behavior far ahead of the rupture, with 



w(x, y) = u(x, y) - \ y y. (47) 

Then 

v 2 ^ = c 2 V 2 w-v(3c 2 V 2 ^ (48) 
ox z ox 

with boundary condition at y = 



dw , ^ d 2 w 
oy oxoy 

w = for x>0. (50) 

Next, Fourier transform along x with 

W(k,y) = J dx exp[ikx]w(x, y). (51) 

Then 
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- k 2 v 2 W = c 2 {^-r- k 2 ){\ + ikv(3)W 



d 2 W 



dy 



fc I 1 — 



dy 2 \ c 2 (l + ikv/3) 
=>. W = W (k)e- 9V , 



W 



g(k) = k< 



c 2 - v 2 + ikvfic 2 ■ i rvA/ \ 

\ — ^ ., 777 with s%) > 

\\ c 2 (l+ikvp) y> 



(52) 

(53) 
(54) 

(55) 



When v < c, in order to insure that the real part of g is positive, one should write it 
as 



g(k) = Vk 2 + e 2 , 



\| c 2 (l + i/b;/3) ' 
where e is small. However, when v > c, one must write instead 



c 2 — f 2 + ikvflc 2 



g{k) = ik t 



v 2 — c 2 — ikvfic 2 
\ c 2 (l + ikv(3) ' 



(56) 



Note that there is a branch of with positive real part everywhere as k moves 
along the real axis. The problem reduces to finding W (k). This function may be 
determined from the boundary conditions. To do so, write 



wq(x) = w(x, y = 0) 



W (k) = f dxw (x)e ikx = [° dxw (x)e ikx = W^{k). (57) 

The superscript — indicates that Wq has no poles in the lower half plane. Next, 
introducing a convergence factor exp[— e\x\] to keep the constant \ y under control, 
sending e to zero at the end of the calculation, write the boundary condition (49) as 



+ V~ e|x| - v(3pf\ e lkx = -gW + + - u/ftA^Wo 

array / e — ike + ik 

-f^S + V- eW -^)^-« + W, (58) 

where the superscript + indicates that Q + has no poles in the upper half plane. 
Therefore, using Eq. (57) one can write 

-g(l + iv(3k)W^ + -\- + = Q + . (59) 

6 %K 6 "T" ^/v 

Define 
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G(k) = g (k)(l + ivPk)/ik 

= ^/(l + iv(3k){v 2 /c 2 - 1-ivkp) = G + (k)G-(k), 



(60) 



with G + (k) = ^v 2 /c 2 - 1 - ivk/3; G~(k) = ^l+iv(3k. (61) 

Note that G + is free of poles or zeroes in the upper half plane (it has a zero in the 
lower half plane) while G~ is free of poles or zeroes in the lower half plane (it has a 
zero in the upper half plane). On the real axis, one takes the branch of both G + and 
G~ that is positive when k = 0; this ensures that the real part of g(k) is positive as 
required. Therefore, write Eq. (59) as 



g+ = -ikG + G~W - + + -A- 

e — ik e + ik 



G+ e-ikG+(0) (e + ik)G+(0) 



ikG~W . (62) 



The two sides of Eq. 62 have poles on opposite sides of the real axis, and must 
separately equal a constant. If the constant is nonzero, then w will be discontinuous 
at the origin, while it if is zero, w is continuous although dw / dx is discontinuous. 
Therefore, take the constant to be zero. One has 



ikW ~ = Xy 1= (63) 

(e + ik)y/l + iv(3k^Jv 2 /c 2 - 1 



\Jv 2 /c 2 



-d 2 w f dk \e~ ikx [°° dk \ y e~ ikx 



dx 2 J 2tt VI + ivfik Jo 2vr jlTwfik 



+ c.c. 



There is a nonzero result only when x < 0. In this case, one must deform the 
contour so that k travels along the positive imaginary axis; k — > ik'. The branch 
of the square root is one that has positive imaginary part on the right side of the 
imaginary axis as one deforms the contour. Making this change of variables, one 
has 

r°° i dk' X v e k ' x 

+ c.c. (64) 



f 

Jo 



lo 2vr VI - vffi 
From to 1/ (3v, the integrand is purely imaginary and cancels with the complex 
conjugate. For the remainder of the contour, the square root is taken on the branch 
with positive imaginary part and gives 



dk' 



A p k ' x 



Ji/v/3 2ir y/v(3k' 
\J-nvf3x 



r°° dk' \ v e k ' x+x/v P 

+ c c =2 / - == 

Jo 2tt V^ 7 



(65) 
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Thus one has 



^V/c 2 - 1 



\ „x/v3 

for x < 



g2w o = j v/^^ 
^ I esse 



dw _ f° ^ X y e x / V/3 
dx Jx ^—irv(3x 



(66) 
(67) 



In particular, the slope of the face of the rupture as x — > — oo is given by 

dwo(— oo) A^ 

ftr ^/c 2 - 1 ' 

In the lab, the slope of the crack face will be 

A,, 



(68) 



\ x Jvy<? - 1 



(69) 



and the opening angle 9 will be 



6 = 2 tan" 1 



A, 



(70) 



A peculiar aspect of the neo-Hookean theory once terms involving E zz have been 
discarded is that it describes material which in its lowest energy state shrinks down 
into a point. This unphysical feature provides an advantage in this calculation, since 
it means that the slope of the crack face is also the same as the slope of the line along 
which material begins to deform as the rupture approaches. As predicted in Section 
2, Eq. (69) is exactly 

In order to determine the velocity of the rupture, one needs a criterion to describe 
when material fails. Consider a bond lying on the center line just before the tip of 
the rupture which in the material frame points along (1/2 \/3/ 2) . In the numerical 
model studied in this paper, the bonds that snap are all of this sort. Experiments 
are carried out in amorphous materials, and it would remain to be shown that this 
type of bond is sufficiently representative of those that snap. Within a continuum 
framework, it is natural to suppose that this bond snaps when 



X 2 



1 2 3 I du 
4 x + 4 I Thj 



(71) 



x=0 
V=0' 



It is easy to come up with more elaborate criteria, but this one has the virtue of sim- 
plicity, and accounts reasonably well both for experimental and numerical results. 
In order to compute the quantity on the right side of Eq. (71), note that 
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/ 



dx 



dw 
dy 



ikx 



-gW = -ik^Wo 



(1 + ikvf3)(e + ik)^v 2 /c 2 - 1 



(72) 



dw 
dy 



dk 



-ikx 



\\Jv 2 /c 2 — 1 — ivkf3 



(73) 



y =o - 271 (1 + ikv(3)(e + ik)^v 2 /c 2 - I 

For x < 0, one must close the contour in the upper half plane, where there are two 
poles, one at i/v(3, and one at ie. These contribute 



dw 
dy 



2ni \ y e x ' v(3 G + (i/v(5) 2m X y G + (0) 



y=0 2tt i v p(-i/ V j3)y/v 2 /c? - 1 2n i^jv 2 /c 2 - 1 



X y e x ^v/c 



-A, 



v 2 1 c 2 



du 
dy 



X y e x ^v/c 



y=0 



v 2 /c 2 



(74) 



and in particular at x = one has 

du 
dy 



1=0 

v=o 



XyV/c 



^v 2 /c 2 - 1 ' 



(75) 



Similarly, for x > one finds that 



<9-u 
dy 



X,,dt 



=0 J 7T 



j-2 e -(t 2 +v 2 /c 2 -l)x/vl3 



(t 2 + v 4)(t 2 + %-i)J^-i 



(76) 



It is also interesting to compute the vertical stress ahead of the rupture, which is 



a y /p& 



= du/dy + f3du/dy 

\ p ~(v 2 /c 2 -l)x/(/3v) 

= K+ y 



^TTX/((3V)^V 2 /C 2 - 1 



Xyjv 2 /C 2 -1 — 



dx > e -(v 2 /c 2 -l)x>/((5v) 



(77) 



Note that while the displacement gradient and strain are finite in front of the rupture, 
the stress does have an inverse square root singularity at the origin. This singularity 
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is due to the Kelvin dissipation, and does not indicate that there is a finite energy 
flux to the tip as in conventional fracture. 

Returning now to Eq. (75), the rupture criterion (71) becomes 



o 1,2 , 3 Ay/c 2 



X f-4 K + 4vyc*-l (78) 

\y = ^ = J I - CyV 2 (79) 



Note that this expression predicts a specific way of assembling samples with dif- 
ferent values of X x and \ y that all should travel at the same speed v. For the exact 
solution of a discrete theory presented in the next section, the final result is of ex- 
actly this same form, with the same quantity appearing, but related to a more 
complicated function ofv/c. 



5.3 Disintegration of back face 



To obtain a final lesson from the continuum solutions, return to Eq. (69). Consider 
two mass points that before the arrival of the rupture lie on the central axis at hori- 
zontal distance dx from each other. According to this expression, on the back face 
of the rupture, they are now separated by the squared distance 

dx 2 + dx 2 X " y 



\i(v 2 /c 2 -iy 

which means that material along the back end of the rupture is stretched by amount 
Aback where 

\2 

\2 __ \ 2 I 

^back ^2 — 1) ' 

Employing Eq. (79), one has that 



A Lk = A*(l - ^) + ~\). (80) 

Inspection of Eq. (80) makes it plausible that extensions along the back end of the 
rupture can be greater than A j : that is, they are generally greater than the extension 
at which rubber near the tip is supposed to give way. This is the reason that material 
must toughen behind the rupture tip. Otherwise, no steady solution is possible and 
the back end of the rupture disintegrates. To emphasize this point, Figure 6 shows 
a numerical rupture solution with bonds in bold when they have stretched beyond 
A/. As predicted by Eq. (80), the entire back surface of the rupture is in this state. 
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Figure 6. Close-up view of top half of rupture with X y = 3.4, X x = 2.2, Xf = 5.5, 
P = 3, traveling atv/c = 1.09. Bonds stretched beyond A/ are drawn much thicker than 
other bonds. As predicted by (80), a line of bonds along the back edge of the rupture has 
been stretched beyond A/, and they remain intact only because of the rule that toughens all 
bonds connected to a node where at least one bond has shrunk below 1.5 times its equilib- 
rium length. Note that several bonds off the main crack line snap as the tip progresses, but 
not enough to destroy the integrity of the material. Without some sort of toughening rule, 
however, the entire back edge of the rupture must disintegrate. 

This calculation explains the need for the toughening rule described after Eq. (33). 
Without such a rule, no supersonic solutions appear in numerical calculations. With 
it they become generic. 



6 Discrete Neo-Hookean theory 

The main weakness in the continuum theory of the previous Section is that the 
rupture criterion is approximate. It is possible to do much better, since one can solve 
analytically the equations of motion for the discrete Neo-Hookean theory given by 
Eq. (37). Recall that this equation was obtained with the following assumptions: 

(1) The coefficient B in (32) vanishes. 

(2) E zz can be set to zero. Because of this assumption, the energy functional is 
quadratic, and the equations of motion are linear. 

(3) Mass points move only vertically. In fact, the horizontal forces on all mass 
points balance, except during a brief time when only one of the bonds has 
snapped for a mass point lying on the crack line. Comparison of analytical so- 
lutions with direct numerical integration of the equations of motion indicates 
that errors introduced by this approximation are on the order of no more than 
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one percent, and a snapshot from a numerical solution of the Neo-Hookean 
theory shown in Figure 7 demonstrates that this approximation is obeyed well 
in the vicinity of the tip. 




Figure 7. Snapshot of supersonic rupture in Neo-Hookean theory taken from numerical 
time evolution of Eq. (37). Note that particles do move almost purely vertically, as shown 
by comparing particle positions with the vertical lines. 

The calculation of steady states for Eq. (37) is lengthy, and relegated to Appendix 
E. The final results are as follows: 

Begin by specifying the dimensionless velocity and damping 

v = v /c, = (3c/a, (81) 

and compute 



c 



3 - cos(cj/£) - 3cj 2 /[4(1 - ifiu)] 
2 cos(ui/2v) 



(82) 



F{u) 



= C + \/C 2 - Iwith abs(0)>l, 

0[iV-l] _ 0-[AT-l] J 

- 2C > cos(a;/2£) 



0^-0- 



-N 



(83) 



and Q(uj) 



F - 1 - cos(cj/2£)' 



(84) 
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Figure 8. Sequence of velocity versus loading curves showing approaches to the continuum 
limit. Curves in (A) correspond to systems of fixed height and fixed continuum dissipation 
(3, but sending the lattice spacing a = L/N to zero. This is achieved through solutions of 
Eq. (85) where N increases and (3 = 0c/ a scales as N. As a goes to zero, the subsonic 
branch of solutions approaches a definite limiting value, but the velocity of supersonic 
solutions increases continually as N increases. Note that on this approach to the continuum 
limit, the fracture energy diminishes as a 2 . Curves in (B) correspond to systems of fixed 
lattice spacing and fixed continuum dissipation 0. As the height L goes to infinity, the 
supersonic solutions approach a limiting value, but the branch of subsonic solutions is 
squeezed into a smaller and smaller region near the origin. 



Then the scaled extension defined in Eq. (79) by 
given as a function of v by 



A,/ V /(4A 2 / - A|)/3 is 



A,, 



vW+I 



exp 




]nQ(u') -lnQ(w') /31n|Q(u/) 



icu'(l + p 2 cu' 2 ) 



+ 



1 + f3 2 uj' 2 
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7 Results 



7. 1 Macroscopic Limit 



Eq. (85) provides a complete expression for the collection of extensions A^ and 
A y that result in a rupture moving at velocity v. Apart from the scaled velocity 
v = v/c, the result depends upon three parameters; the system height N, the exten- 
sion A/ at which bonds snap, and the coefficient of Kelvin dissipation (3 = c(3/a. 
One can now search for the conditions under which one should expect subsonic 
fractures, and the conditions under which one should expect supersonic ruptures. 
First, consider systems of fixed height and fixed level of dissipation (3 as the lattice 
spacing tends to zero. This situation is described by fixing L and taking the limit as 
N — > oo of Eq. (85) with (3 = c(3N/ L also scaling as N. Figure 8 (A) shows that in 
this limit, there is a narrow band of subsonic solutions followed by a broad band of 
supersonic solutions. Note that since the failure extension A/ remains fixed while 
the lattice spacing a = L/N vanishes, the fracture energy vanishes as a 2 during 
this limiting procedure. Thus, this limiting procedure, which at first seems the most 
sensible, corresponds to something physically rather odd. Alternatively, one can set 
f3 to a constant and send N to oo so that the sample becomes infinitely high. In this 
limit, plotting solutions versus X y , the subsonic ruptures disappear, and only super- 
sonic solutions survive, as shown in Fig. 8 (B). However, ones conclusions about 
the true nature of this macroscopic limit depend upon how one scales the solutions, 
as illustrated in Fig. 9. For a system of any given height, there are both subsonic and 
supersonic solutions. The subsonic solutions are found at small strain, and as N be- 
comes large, the range of extensions \ y that produces them becomes progressively 
smaller. There is a plateau near the wave speed that becomes wider and wider as 
N increases. Finally, supersonic ruptures appear for extensions \ y on the order of 
\f. The point to emphasize is that depending how extensions are scaled, either the 
supersonic or subsonic branches can be viewed as the macroscopic limit. In most 
brittle materials it is impossible for cracks to reach the wave s peed because they be- 
come unstable to side-branching before this point is reached dFineberg and Mardei , 
1999). One of the things that appears to make rubber different is that the ruptures 
are so stable that it is possible for them to pass the wave speed and move beyond it 
without instabilities intervening. 



7.2 Comparison with Experiment 



To close this investigation, I compare the results with experiments on rupture of 
rubber sheets. It was already demonstrated in Figure 3.2 that the Mooney-Rivlin 
theory adequately captures the varia tion of sound s peed with extension. The re- 
maining two quantities measured by Pet ersan et al.1 d2004) are rupture speeds and 
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Figure 9. Four different views of Neo-Hookean crack velocities, showing that depending 
upon how they are scaled and displayed, one focuses either upon conventional subsonic 
fractures, or supersonic ruptures. The definition of X y is given in Eq. (79). In the limit of 
infinite system height N, the two different types of solutions are separated by an infinitely 
long plateau at the wave speed. 



opening angles. Before making the comparison, two reasons to view the compari- 
son with a bit of skepticism should be noted. First, rubber is an entangled polymer 
network, not a triangular lattice. Second, although the Kelvin dissipation propor- 
tional to (3 plays a very important role in the theory, no estimate of its value from 
experiment has been provided. The reason is that dissipation in real amorphous 
solids is not of the form employed here, nor is there any simple way to correct the 
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deficiency. The spatial decay rate of sound waves in rubber over freq uencies ran g- 
ing from kilohertz to megahertz is almost perfectly linear in u (iMott et all 12002): 



a = Cm, 



(86) 



where C is a constant. However, given the form of Kelvin dissipation employed in 
this paper, sound decays at the rate 



It is impossible to find a value of (3 that makes Eq. (87) a good fit to Eq. (86). 
A more realistic rule for Kelvin dissipation would provide a frequency-dependent 
sound speed according to 



the form of dissipation used in Eq. (41) corresponds to sending to the high-frequency 
sound speed to infinity. However, the frequency dependence of sound attenu- 
ation does not resemble experiment any better after inclusion of c^. Thus I will 
simply contin ue to use the s i mple st form of Kelvin dissipation, as it is familiar and 
conventional dFradkin et all 120031) and take the dimensionless measure of dissipa- 
tion, (3 = (3c/ a to be of order unity. Fortunately, none of the final results depend 
much on the value of (3. 

Figure 10 assembles experimental and theoretical results. According to the theory 
for triangular lattices, samples with extensions X x and \ y depend only upon the 
scaled variable \ y given by Eq. (79). This scaling of the velocity is compatible 
with all the data. Thirteen experimental trials where ruptures ran straight collapse 
onto five points, with rather little variation in the scaled velocity. The scatter in 
the data is rather large, and therefore consistent both with the simplified results of 
Section 5, as well as the more elaborate results of Section 6. The figure also shows 
a comparison of direct integration of the equations of motion, Eq. (33). For the 
equations of motion in this figure, the Mooney-Rivlin parameter B has been set to 
zero, but E zz has not been eliminated from Eq. (17). Agreement with the analytical 
results from Eq. (85) is excellent, showing that details of how rubber relaxes behind 
the tip of the rupture do not have much effect on the dynamics. As already shown in 
Figure 5, rupture speeds are not measurably affected by including B. Therefore, the 
analytical results of Eq. (85) capture rupture speeds, experimental and numerical, 
rather completely. The simple result of Eq. (79) is adequate for a first pass. In fact, 
the value A/ = 5.5 of the bond failure extension was obtained by fitting (79) to the 
experimental data, and this value of A/ was then used unchanged in all numerical 
runs. 

Finally, Figure 11 compares experimental results for crack opening angles with 
predictions based upon numerical solutions of the most realistic numerical system, 




(87) 




(88) 
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(33). That is, the nonlinear terms from E zz that appear as rubber shrinks towards its 
equilibrium relaxed state, and Rivlin's nonlinear contribution to the Mooney-Rivlin 
energy are all included. Agreement between theory and experiment for the open- 
ing angle is still not completely satisfactory. The experimental points are widely 
scattered, indicating that the reduction to the variable X y may not be appropriate, 
and experimental values lie systematically below theoretical predictions. Either the 
simplistic form of the dissipation, or the simplistic triangular micro structure might 
be to blame for this discrepancy. 



There is one final potential discrepancy with experiment that should be mentioned. 
According to the theory, the dynamical solutions do in clude subsonic ru ptures at 
small extensions. Many rubbers are well known to creep (IHui et al.L l2003). and tear 
slowly in trouser tests, but in our biaxially loaded samples of natural latex rubber 
we never observed cracks to creep, or to travel slower than the sound speed at all. 



8 Conclusions 



The main points established by this theory for the rupture of rubber are the follow- 
ing: 

(1) The rupture of rubber is a shock phenomenon, with the back edges traveling 
at a wave speed, and the tip of the rupture consisting in the place where two 
shocks meet at a point. 

(2) The essential physical ingredients are dissipation and some toughening that 
allows the back end of the rupture to retain its integrity. Static hyperelasticity 
appears not to be relevant. 

(3) Predictions for rupture velocity are in satisfactory agreement with experiment. 
Predictions for opening angle are less so, perhaps because the computations 
have been performed in a triangular lattice, and with a very simple form of 
dissipation. 

Additional physical question that have not yet been resolved are 

(1) Under what conditions do cracks in rubber creep, and when instead are they 
supersonic? 



(2) What is the origin of rupture path oscillations reported by Dee gan et al.l (12002)? 
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Figure 1 1 . Rupture opening angles obtained from simulations and compared with experi- 
ment. The agreement is not very satisfactory; the experimental results are widely scattered, 
but lie systematically below the numerical predictions. For prediction of this quantity, it 
may be that treating rubber as a triangular lattice is not adequate. 



A Reduction to 2 dimensions 



This Appendix shows in two different ways how to obtain an effective two-dimensional 
equation of motion for a rubber sheet. In the first, method, the incompressibility of 
rubber is used to calculate the thickness of the sheet at every point, thus express- 
ing displacements across the thickness of the sheet (z direction) in terms of the 
extensions along x and y (Figure 1). 



A.l First method 



Rubber is highly incompressible, so one can set 



det 



du 



df 



0. 



(A.1) 



For a thin sheet, assume that one can neglect du y jdz and du x /dz, on the grounds 
that u x and u y should be uniform through the thickness of the sheet. Then Eq. (A.l) 
becomes 

du z ( du x du y du x du y \ 



dr z \dr x dr y 



dr y dr x 



Since u x and u v are assumed to be independent of z, one can write 

-l 



du x du y du x du y 



■ur 



dr x dr y dr y dr x 



(A.2) 



which is odd in r z In moving to a two-dimensional theory, replace all quantities by 
their averages across the sheet. That is, if the sheet has thickness t in the reference 
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frame, then for example, 



E xz (r x , fy) — 



*/2 dr z 

-t/2 t 



(A.3) 



Consider now 



E xz {r) 



1 / du x du x du y du y du z du 



+ 



+ 



2 \ dr x dr z dr x dr z dr x dr 



To obtain the two-dimensional version of this quantity, note that the first two terms 
vanish because of the derivatives with respect to r z while the last term is odd in r z , 
and vanishes when averaged across the sheet thickness in (A.3) . Therefore, in the 
two-dimensional theory, one can take E xz = E yz = 0. Finally, consider 



Ely? = — 



du" 
dr, 




du x du y du x du y 
dr x dr 1t 



-2 



dr y dr x 



- 1 



(A.4) 



The derivatives appearing in the denominator of (A.4) can be expressed in terms of 
the two-dimensional invariants in Eq. (7) as 



Ezz 2 V [/, -i 2/, 4 1 

which is Eq. (8). 



A.2 Second Method 



Specialize to the case of a Neo-Hookean material. Note that raised indices are em- 
ployed on u because elsewhere in the manuscript subscripts are needed to index 
the locations of multiple particles; raised indices are just ordinary Cartesian com- 
ponents of the vector u. For an incompressible solid the Cauchy stress tensor is 
( Ogden[fl984l) 

o du a du^ 



T, 



a(3 = PC 



<9r 7 dr 1 



-p6 a 



Pi 



(A.5) 



where p is a pressure that must be determined by the condition of incompressibility. 
To find an equation of motion, one needs the Piola-Kirkhoff stress tensor, which for 
an incompressible material takes the form 



Sat 



dr a 
du x 



T: 



A/3 = pc 



, du 13 dr a 



P 



(A.6) 



Writing out the equation of motion gives 



d 2 u c 



_d_ 

dry 



S. 



Xa 



_d_ 

dr x 



pc 



,du c 
dry 



dr x 
du a 



(A.7) 
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From Eq. (A. 2) one can write dr a /du^ as 



/ dr^_ dr y 

du x du x 

dry_ 

\ du y duv 



du x du v 
dr x dr x 

du x du y 



-1 



du z 
dr. 



duy 

dr v 



dr x 



du x du x 
dr y dr x 



(A.8) 



Now examine the equation of motion for u 1 



d 2 u x 
W 



P- 



d 2 du a 
dr x PC dr x 



d ( du z du y \ d ( du z du y " 



dr x \ dr z dr. 



dr y \ dr z dr a 



(A.9) 



Compare this result with the one that would come by inserting the constraint at the 
outset: Using Eq. (A. 4) one finds 



d 2 u x 
W 



P- 



Su x (r) 



- / dfpc 2 (E xx + E yy + E z 
A J 



' _d_ 

dry 



-pc 



,du c 
dry 



_d_ 

dr x 

_d_ 

dr„ 



pc 



pc 



' du z \ du y 



\ dr z J dr y 
du z \ 3 du y 



(A. 10) 



dr 7 } dr n 



Eqs. (A.9) and (A. 10) are the same provided that 



p = pc 



' du z 
dr 7 . 



(A.ll) 



Thus the equation of motion obtained by employing the constraint in Eq. (8) is 
compatible with the equation of motion one obtains from the Piola-Kirkhoff stress 
tensor so long as one uses Eq. (A. 1 1) for the pressure. Furthermore, this expression 
for the pressure is precisely what is needed so that T zz vanishes in Eq. (A. 5), and 
that in turn is what one would expect as the appropriate boundary condition for a 
thin sheet. 



B Sound Speeds 



Given an energy functional 



U = pj dre(h,I 2 ), 
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where p is the mass per area measured in the reference frame, the aim of this ap- 
pendix is to calculate soun d speeds. The same results are found for example in 
(Eringen and Suhubl 19741 v. 1, pp. 120, 263), but to obtain the simple expres- 
sions for longitudinal and shear waves needed here, it may be easier to begin again 
than to work backwards through so much notation. To begin, find how U varies 
when there is a small change in u: 



1 8U ^ 



d 2 u 1 de 1 
dr a drp dE a p 2 



du 1 d du 1 d 



dr a dr 



+ 



de 



P 



drp dr a J BE, 



af3_ 



Take de/dE a p to be symmetric under interchange of a and (3. This is only true if 
from now on whenever one sees E xy in some term in the free energy, one replaces 
it by (E xy + E yx )/2, and we will have to be careful to do that. However, assuming 
this symmetry, one can write 



E 

a,P 



1 su 

p8ift{r) 
d 2 u< de 



E 

a,/3 



d 2 u< de 
dr a drp dE a 



dr a drp dE, 



ap 



^ 2 dr 



-E 

a,/3 



d 2 vF 



xP 

'A 

drp 
de 



+ E 



du* dE a < 



P' 



d 2 e 



a 7^, dr a drp dE a pdE a ,p, 



du 1 du 1 
dr a i drpi 



d 2 e 



dE a pdE a ipi 



dr a drn dE a 



xP 



+ E 



du 1 du 1 ' <9 2 ti 7 ' 



(B.l) 



d 2 e 



Jpfy dr a dr a i drpSrp dE a pdE a ,p> 



Now take u to represent a sheet loaded up in biaxial strain, and superpose a small 
amplitude wave with polarization e traveling with wave vector k. Take A 7 to be the 
extension factor along direction 7, and r 7 to be a position coordinate in the material 
frame. Keep only terms of order e. Inserting such a plane wave into Eq. (B.l), the 
result is 



1 SU 



P6u-r(r) tp 

=E 

a,P 



a ' kp> kpey — — — — 

OE a p q//3 , 7 , OEapOEa'P 1 



k a kpt-y 



de 



dE, 



: E 

/3f3'Y L 



kpkp>e 1 



af3_ 



de 



*S X-y X^kpikpe^ 



d 2 e 



x de 



+ Xj XykpikpSj/ 
d 2 e 



dE 1 pdE 1 ipi 
d 2 e 



dEpp dE^pdEypi 



7 dEpp 



^ ^ dEypdE^i pi 



kp'kpCy. 



(B.2) 
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Therefore, one has an equation of motion for sound waves 



poj 2 e 1 



5U 
8u' y (r) 

j>W L 



Be 



BE, 



d 2 e 



(B.3) 



B. 1 Specific expressions for longitudinal and shear waves 



Take k = k(l, 0). Assume that 

d 2 e d 2 e 



BE BE BE BE 

^■^xx^-^xy ^-^yy^-^xy 



0. (B.4) 



This will always be the case in biaxial strain just so long as the energy only depends 
upon strain through the combinations in 1\ and I 2 . Longitudinal waves are found 
by looking for a wave polarized along x, which means that only e x is nonzero, so 
7 = 7' = x. Note in addition that in Eq. (B.3) one can have kp — kp> — 1 only if 
(3 = (3' = x. Therefore 

, Be x9 B 2 e ^ 

Next look for a wave polarized along y. Now 7 = 7' = y. One still has to have 
(3 = (3' = x. Therefore for shear waves 



Be B 2 e 

c 2 x = 7^ + >? y -QgT Assuming symmetry (B.6) 

It is easy to use Eq. B.6 improperly. It is only valid if e is treated as a symmetrical 
function of E xy and E yx , and if partial derivatives with respect to these two quanti- 
ties are independent. It is hard to remember to retain this convention, and it is safer 
simply to set E xy = E yx and treat e just as a function of one of them. In this case, 
one must write 



n Be \ 2 V B 2 e m „ N 

The analogous expressions for wave speeds along y follow by flipping the roles of 
x and y. 



The expressions for sound speeds take much simpler forms if one establishes spa- 
tially uniform states in the rubber and considers small uniform disturbances. Return 
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to Eqs. (15a) and (15b). Establish the displacement field 



u = (X x x + s xy y, X y y + s yx x) . (B.8) 

Specialize now to the case where a sample is subject to uniform bi-axial strain, and 
sheared in the y direction, so that s xy is zero. Then inserting Eq. (B.8) into Eq. 5 
gives 



E x 
E, 



yy 



2 {^x + s yx l) 



E =E 

J -'xy J -'yx 



IX B 

2 y y x 



(B.9) 
(B.10) 
(B.ll) 



Derivatives with respect to components of the strain tensor can all now be expressed 
in terms of the new variables \ x , X y , and s yx . Since these variables correspond 
exactly to quantities one controls experimentally, it is good to express sound speeds 
in terms of them. One has 



d \E xx E yy E xy } 
d {X x X y s yx } 



<K o O 

Xy 



^ s yx /2 X y /2 



(B.12) 



Inverting this matrix, one has 



d {X x X y s yx } 



xx&yy&xy} 



Therefore, one can write 



( i_ 




s 2 
yx 

X X Xy 



0_ v x 

Ay 



^Syx 
X x Xy 



2_ 

Xy 



(B.13) 



d 1 d 

dE xx X x d\ x 



(B.14) 



d 2s vx d 2 d m _ 

- + t — — • (B.15) 



dE xy X y X y dX x X y ds yx 



Inserting Eq. (B.14) into Eq. (15a) and Eq. (B.15) into Eq. (15b) and evaluating at 
= gives Eqs. (14). 
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C Force computation 



I record here some methods used to calculate forces in molecular dynamics that 
assist in writing computer code, and that may not have been published previously. 
The force on component a of particle I is defined to be 

dU 

s -esp (C1) 

For the purposes of writing computer code, it is not efficient to proceed directly with 
this expression, because in the course of computing, say, F\ one might calculate 
some quantities that will also appear in F 2 and efficient code will duplicate as little 
computation as possible. Therefore, define the insertion operator Any term that 
is multiplied by X; is to be inserted into the memory location that holds the force 
on particle /. So one computes 



dw dl{ dw dL 



+ 



dl\ duf dl\ duj 



(C.2) 



Thus one must compute a sum of two terms. The first is 



dw 1 

li,jen(i) 1 i,jen(i) 



(but this term vanishes if u^ > aXf) 
The second is 

dw 3 



E te-z,- 



dw 1 
dl{ 3a 2 



u 



E fr-^K 

ijen(i) 



13 dJh 4a 4 



2 4 

-Ft - -{uij ■ - a 2 ) 



(C.3) 



(C4) 



(but vanish if u^ > a A/) or, if one takes the alternate representation of I\, 



E Pi-ZiK 



dw 9 



-{Uij -Uij -a ) 



+ E3 



dw 9 dHj 
dl\ 8a 4 duf ' 



(C.5) 



The final term to compute is (with g ijk = (m^ • u ik + 2a 2 ) 2 ^ 



1 1 

H * = 27 51 ^ ' " ifc + 2a2 YK u ij)K u ik) = ^ E 9ijkhijh ik (C.6) 



j^k£n(i) 
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and derivatives of this object contribute to the force 



- E 

i 



dw 1 9 



h'ijhikgijk^- (Xj — Xj) + ^ik^-ijdijk^fr ~ Xfc 



I'. J'.hU;,!,- (?i {ufj + < fe } - XfcUg - X 



> . 



(C.7) 

Note in all these expression that various quantities need only be computed once and 
inserted into registers for particles i, j, and k. The insertion operators X , lj , and 
l k keep track of which quantities to put where. 



D Lattice Instabilities 



I found numerically that when I employed Eq. (28), the uniformly strained lattice 
would spontaneously develop a striped pattern when stretched beyond a critical 
value . To analyze this problem, I computed the phonon dynamical matrix (M ardei . 
2000, Eq. 13.8, p. 307). The calculation is a straightforward exercise in phonon 
physics, and no details need to be reported. Negative eigenvalues of this matrix 
indicate instability of the uniform state. As shown in Figure D. 1 , for uniform biaxial 
strain a bit above \ x = \ y = 3, the spatially uniform lattice becomes unstable. 
However, if one employs instead Eq. (29), then as shown in Fig. D.2, the lattice 
remains stable. For this reason, Eq. (29) was usually employed, despite its greater 
numerical cost. It would be interesting to see whether the instability in Fig. D. 1 is 
related to strain crystallization. 
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Figure D.l. Frequency squared as a function of wave number for various levels of uniform 
biaxial strain, using the energy potential function in Eq. (28). As the strain increases, the 
lattice becomes unstable to a distortion in which every other column of atoms moves in 
opposite directions, similar to motions of atoms in optical modes, or the Peierls distortion. 
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Figure D.2. Frequency squared as a function of wave number for various levels of uniform 
biaxial strain, using the energy potential function in Eq. (29). With this representation of 
the strain invariant, the uniform lattice remains stable against high-frequency distortions 
over the range of extensions employed experimentally. 
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E Solution of Discrete System 



This Appendix contains the steady-state solution of 



jen(«) 



The methods employed are tho se o f|Slepvanl(fl98llll982Ll2002h : see also (Ma rder and Grossl 



1995t lMardeiil2004l) . ISlepyanl(l2002L p. 4781 notes the existence of supersonic solu- 
tions for a related problem. However, the steps below are worth recording in detail 
because the particular combination of Kelvin viscosity, Mode III, and a strip of fi- 
nite height needed here has not been published, although no especially new ideas 
are involved. 

Move to notation that explicitly describes locations in a triangular lattice, u(m, n), 
where u describes the vertical motion of atoms only, since the horizontal motion 
is neglected, and where m 6 (—00 ■ • • — 3, —2, —1, 0, 1, 2, 3, . . . 00) and n G 

( 5/2, -3/2, -1/2, 1/2, 3/2, 5/2 ... . Then in steady state, one has the 

symmetries 



u(m, n, t) —u(m+ l,n, t + a/v) (E.2) 
u(m, n, t) = —u(m, —n, t — a[l/2 — g n ]/v) (E.3) 
u(m, 1/2, t) = -u(m, -1 /2, t - a/2v), (E.4) 



where 



1 « „ = 3/2, 7/2... 

' ^0 ,ifn= 1/2, 5/2... 



Assuming that a crack is in steady state, we can therefore eliminate the variable m 
entirely from the equation of motion, by defining 

u n (t) = u(0,n,t). (E.6) 

Next, define dimensionless variables 

t — tc/a; (3 = (3c/ a, and v — v/c. (E.7) 

Then 
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tin(f>f(l + ^) 



+u n+1 (t- (g n+1 - l)/v) +u n+1 (t - g n+1 /v) 

+u n (t+l/v) -Qu n (t) +u n (t - l/v) 

+tt n _l(f- (# n _i - l)/u) +Un-l(* - 9n-l/v) 



if n > 1/2 and 



tii/2 (*) = 



+(i + / g|) M3/2 (f) + (i + ^) U3/2 (t- i/v) 

+ (1 + /3§) [u 1/2 (t + l/v) - 4« 1/2 (i) + «i /2 (f- 1/5) 

+fl(-*)(l+^fl)[«_ 1/2 (f)-u 1/2 (t)] 

+0(1/(26) - t)(l + /3f )[u_ 1/2 (f- l/v) - u 1/2 (t)\ 



(E.8) 



ifn = 1/2. 



The time at which the bond between u(0, 1/2, t) and u(0 : — 1/2, t) breaks has been 
chosen to be t = 0, so that by symmetry the time the bond between u(0, 1/2, t) and 
u(l,-l/2,i) breaks is 1/26. 

Above the crack line, the equations of motion are completely linear, so it is simple 
to find the motion of every atom with n > 1 /2 in terms of the behavior of an atom 
with n — 1/2. Fourier transforming in time gives 



— uj 2 u n {uj) = -(1 — if3ui) 



U n+1 {uj) [ e M9n+l-l)/v + giw^n+ij/ti] 

+ u n (to) [e iuJ ^ _ 6 + e -^/ c ] 

+ Un-!(u) [e^n-i-lj/v + giwO/n-O/ti] 



(E.9) 



Let 



3 fc(fi-l/2)-iwfl„/(25) 



u^wj = u 1/2 [u)e^" ~<-' ~»»'v-/. (E.10) 
Substituting Eq. (E.10) into Eq. (E.9), and noticing that g n + g n+1 = 1 gives 



-u; 2 -ui/ 2 (cu) = (l—i/3u)- 



u 1 / 2 {uj)e k ^ e iu} (9n+i+g n -2)/Civ) _|_ e iu>(g n+1 +g n )/(2v)^ 

+ U1/2H [ e W« _ 6 + e"^/ 5 ] 

+ ■u 1 / 2 (u;)e~ fc [e iw ( s "- 1+s "- 2 )/( 2S ) + e iw ( s "- 1+s,n )/( 2C )] 

(E.ll) 



cj 2 4 

+ - [2 cosh(fc) cos(w/(2i5)) + cos(u/v) - 3] = 0. (E.12) 



l-i/3u 3 
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Defining 



3-cos(a;/i))-3a; 2 /[4(l-^)] 
2 cos(cj/2t>J 



one has equivalently that 



= e fc = C + ^C 2 - 1 with abs(0) > 1. (E.14) 
One can construct a solution which meets all the boundary conditions by writing 

jAN+l/2-n] _ .-[tf+l/2-n] 

^iM-K^F ] 

iv + in 1 

This solution equals ui/ 2 for n = 1/2, and equals XJ^ej (e 2 + uj 2 ) for n = N+l/2. 
The reason to introduce e is that for n = N + 1/2, w(m, n, t) — Un- The Fourier 
transform of this boundary condition is a delta function, and hard to work with 
formally. To resolve uncertainties, it is better to use instead the boundary condition 

u N+1/2 (t)=U N e-^, (E.16) 
and send e to zero the end of the calculation. In what follows, frequent use will be 
made of the fact that e is small. 

The most interesting variable is not uy 2 , but the distance between the bonds which 
will actually snap. Furthermore, the quantity multiplied by the 9 function in Eq. 
E.8 is operated upon by (1 + f3d/dt) since dissipation stops operating when bonds 
break. For this reason define 

u ^ = u 1/2 (t) - u_ 1/2 (t) = u 1/2 (t)+u 1/2 (t + l/2v) (E n) 

W(t) = (l+^)U(t) 

Rewrite E.8 as 



U l/2(t) = - 



+(1 + (3i)u 3/2 (t) + (1 + (3f t )u 3/2 (t - 1/v) 
+(l + /3f) (u 1/2 (t + l/v)-A + u l/2 (t)+u 1/2 (t- 1/v)) 
-2U(t)9(-t) - 2U(t- l/2v)9(l/(2v) - t) 



(E.18) 

Fourier transforming this expression using Eq. (E.15) and defining 

U ± (u)= dte tuJt U{t)9{±t), (E.19) 

J — oo 
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W ± (u) = / die iuji W(t)9(±i) 

now gives 

(1 - iup)u 1/2 (u)F(u) - (1 + e iuj/2 *)U-(uj) 

with 



U N 2e 
~N~ oo 2 + e 2 ' 



F(u) = { — — — 2Qcos(u;/25) + 1 



Next, use Eq. (E.17)in the form 



W(u) = (1 - ^(w) 



to obtain 

Writing 
finally gives 

with 



W(u)F(u) - 2(cos 2 u/4v)U-(uj) 



U N 2e 
~/V~ u 2 + e 2 ' 



W(u) = W + (u) + W~(u) 



W + {u)Q{u) + W~{u) = U N Q 



1 1 

+ 



e + iu e — iu 



(E.20) 

(E.21) 
(E.22) 

(E.23) 

(E.24) 
(E.25) 

(E.26) 



Q = F/(F-l-coa(u/2v)). 
The Wiener-Hopf technique dNoblelll958h directs one to write 



Q(u) 



Q+(u) 



(E.27) 



(E.28) 



where Q~ is free of poles and zeroes in the lower complex u plane and Q + is free of 
poles and zeroes in the upper complex plane. One can carry out this decomposition 
with the explicit formula 



Q ± (u) = exp[\imJ 



du' lnQ(u') 
2ir iu =)= e — iu' 



(E.29) 



Separate Eq. (E.26) into two pieces, one of which has poles only in the lower half 
plane, and one of which has poles only in the upper half plane: 



W + (u) Q U 



N 



1 



QqU 



N 



W~(u) 



Q+(u) Q(0)(-iu + e) Q-(0)(iu + e) Q~(u) 



(E.30) 
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Because the right and left hand sides of this equation have poles in opposite sections 
of the complex plane, they must separately equal a constant, C. The constant must 
vanish, or U~ and U + will behave as a delta function near t = 0. So 



\*r-( \ tt ( w ) a ixr+< \ tt QoQ fa) 

W (u) = U N , and W + (u) = U N — . (E.31) 

Q (0)(e + iuj) Q {0){e — iu) 

One now has an explicit solution for W(u). Numerical evaluation of W(t) from Eq. 
(E.31) is fairly straightforward, using fast Fourier transforms. The most interesting 
quantity to obtain is the separation between bonds opposite the crack line at t — 
0, since by setting this quantity so that the bond snaps, one obtains a consistent 
equation of motion. So one wants to find Z7 (0). To obtain it, write 

= njiEM. (E32) 

1 — ij3uj 

The denominator of Eq. (E.32) has a pole in the lower half plane at 

uj — —%j /3 = — iujq (E.33) 
and this pole must be subtracted off in order to form U~ . So one has 

= H^M-ty^o) (E34) 

1 — i(3u! 

In order to find U(t — 0) it is sufficient to find 

U(t = 0) = lim iuU~(u;). 

The reason is that U~(t) is zero for all positive t, dropping to zero right at t = 0. 
The value of U(t = 0) is given by the discontinuity in U~(t). However, if U~(uj) 
decays faster than l/uo for large lu, then the inverse Fourier transform of U~{u) 
must be continuous at t = 0. Therefore, from the coefficient of l/iu, one can pick 
out the value of U(t = 0). Since W~(t) like U~(i) has a step-function discon- 
tinuity at t — 0, W~(uo) decays as 1/iuo as u goes to infinity. Thus one deduces 
immediately from (E.34) that 

U(t = 0) = co W-(-iuj ). (E.35) 

Returning to (E.31) one has 



TT TTU H\ TT (— ^0) ~~ 

Uo = U{t = 0) = C/jv • (E.36) 

Note that the height of the top of the system [/ ^ obeys 

^-a(N + l/2)\ y = U N , 
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and that Q = Q(0) = 1/(2N + 1) so 




The condition for a bond to snap is that the total length of the bond reach aXf. Note 
from the definition in Eq. (E.17) that U gives only half the bond length. Therefore 



Employing (E.29) to obtain an explicit representation of Q leads to Eq. (85). 
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